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We derive an infinite liierarcliy of exact equations for tlie Bak-Sneppen model in arbitrary dimen- 
sions. Tliese equations relate different moments of temporal duration and spatial size of avalanches. 
We prove that the exponents of the BS model are the same above and below the critical point and 
express the universal amplitude ratio of the avalanche spatial size in terms of the critical exponents. 
The equations uniquely determine the shape of the scaling function of the avalanche distribution. 
It is suggested that in the BS model there is only one independent critical exponent. 
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Recently Bak and Sneppen ^ introduced a particu- 
larly simple toy model of biological evolution (BS model). 
The model correctly reproduces such features of real evo- 
lution process as punctuated equilibria, power law proba- 
bility distributions of lifetimes of species and of the sizes 
of extinction events. In spite of the simpHcity of the 
rules, the model exhibits extremely rich and interesting 
behavior. 

In the one-dimensional version of the model L numbers 
fi are arranged on a line. At every time step the small- 
est number in the system and its nearest neighbors are 
replaced with new uncorrelated random numbers drawn 
from the uniform distribution between and 1 . The gen- 
eralization to higher dimensions is straightforward. In 
fact there exists a whole class of models where one selects 
the site with the extremal (global maximal or minimal) 
value of some variable and then changes this variable 
at this site and its nearest neighbors according to some 
stochastic rule. One of the best known representatives 
of this class is Invasion Percolation ||^. Such models, re- 
ferred to as extremal models were extensively studied (for 
a recent review see |^). 

The interesting feature of the BS model (as well as 
other extremal models) is its ability to organize itself into 
a scale-free stationary state. The dynamics in this critical 
state is given in terms of bursts of activity or avalanches, 
which form a hierarchical structure of subavalanches 
within bigger avalanches. Here we introduce a "master" 
equation for this avalanche hierarchy. It describes the 
cascade process of smaller avalanches merging into big- 
ger ones as the critical parameter is changed. From this 
equation we derive an infinite series of exact equations, 
relating different moments of temporal duration S and 
spatial size of individual avalanches. 

The "master" equation connects undercritical and 
overcritical regions of parameters. Given the existence 
of the scaling, we rigorously prove that the exponents of 
the BS model are the same above and below the transi- 
tion. From our results it follows that all terms of Taylor 



series of the scaling function f{x) for the avalanche dis- 
tribution are uniquely and explicitly determined by two 
critical exponents of the model. It was suggested that 
usual restrictions on the shape of f{x) indirectly relate 
these two exponents and, therefore, reduce the number 
of independent critical exponents in the BS model to just 
one. 

To define these avalanches one records the signal of 
the model, i.e. the value of the global minimal number 
fmin(s) as a function of time s. Then for every value of 
the auxiliary parameter /o, an /o-avalanche of size (tem- 
poral duration) S is defined as a sequence of S* — 1 suc- 
cessive events when f minis) < fo confined between two 
events when f minis) > fo- In other words, the events 
when f minis) > fo divide the time axis into a series of 
avalanches, following one another. It is easy to see that 
an avalanche defined by this rule is nothing else but a 
creation-annihilation branching process where sites with 
fi < fo play the role of particles. The avalanche is ter- 
minated (and the next one is immediately started) when 
there are no such "particles" left in the system. As in any 
other creation-annihilation branching processes (such as 
directed percolation) in BS model there exists a criti- 
cal value fc of the creation probability fo, for which the 
creation of particles is just marginally balanced by their 
annihilation, and avalanches of all sizes can be realized. 
In the stationary state of the BS model on the infinite 
lattice, f minis) < fc for every s. Therefore, the overcrit- 
ical ifo > fc) region of the branching process parameters 
is not accessible, since there are no events starting or 
terminating such avalanches. However, if the system is 
artificially prepared in the overcritical state with fi > fo 
everywhere, one can observe overcritical avalanches. In 
this case there is a non-zero probability Poo(/o) to start 
an infinite avalanche and at the same time the size of 
finite avalanches has a finite cutoff. 

The events within the same avalanche are spatially and 
causally connected. It is easy to understand that the po- 
sition of active site at any time step within the avalanche 
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is connected to the set of sites covered (updated at least 
once) by the avalanche up to this time step. We charac- 
terize an avalanche by two principal numbers: 1)5'- the 
avalanche size, equal to its temporal duration; 2) ricov ~ 
the number of covered sites, i.e sites that had their ran- 
dom number updated at least once during the course of 
this avalanche. In one-dimensional models the connected 
nature of the set of covered sites ensure its compactness 
and, therefore, ricov is equal to the avalanche spatial ex- 
tent R. In higher dimensions (below the upper critical 
dimension) it was conjectured in |^ that the set of cov- 
ered sites is a non-fractal object of the same dimension- 
ality d as the underlying lattice. In this case the spatial 
size R of the avalanche can be defined by the relation 
n —7?'^ 

"'cov — -'I' ■ 

/o-avalanches in the Bak-Sneppen model were shown 
to be exactly equivalent to the realizations of the BS 
branching process Q . In this process one only keeps track 
of the numbers fi < fo, and at each time step activates 
the smallest one of them. The BS branching process is 
terminated when there are no numbers fi < fo left in the 
system. Besides the fact that the BS branching process 
is a very effective way to simulate the BS model numeri- 
cally, it has the additional important advantage that the 
overcritical region fo > fc becomes accessible. 

The quantity of primary interest in the BS model is the 
probability distribution P{S, fo) of the avalanche sizes S 
at any given value of the auxiliary parameter fo- The 
moments in time, when fo < f minis) < fo + dfo serve 
as breaking points for /o-avalanches but not for fo + dfo- 
avalanches. Therefore, when is raised by an infinitesi- 
mal amount dfo some of /o-avalanches merge together to 
form bigger {fo + (i/o)-avalanches. In the rest of this pa- 
per we study in more detail the properties of this merging 
process and the avalanche hierarchy that it induces. 

The most important observation about /o-avalanches 
in the BS model (as well as in several other extremal 
models, such as the Sneppen model or Invasion Per- 
colation Ql) is that when an /o-avalanche is terminated, 
the numbers fi on the set of Ucov — R'^ updated sites are 
uncorrelated and uniformly distributed between fo and 1. 
To prove this statement one notices that at any single 
time step all new random numbers are drawn from the 
uniform distribution between and 1. If such a num- 
ber happens to be larger than /o, its exact value is not 
important for the dynamics of this avalanche. 

The direct consequence of this observation is that the 
probability of an /o-avalanche of spatial size R'^ to merge 
with the subsequent one when the parameter fo is raised 
by dfo is given by R'^j^^. (the merging occurs if at 
least one of the changed numbers falls in [fo, fo + dfo]-) 
For the following arguments to be true it is important 
that any two subsequent avalanches are mutually un- 
correlated. That is: the probability distribution of /q- 
avalanches, starting immediately after the termination 



of an /o-avalanche of a given size S is independent of S. 
That is true for the BS model since the dynamics within 
an /o-avalanche in BS model is completely independent 
of the particular value of the numbers fi > fo in the 
background that were left by the previous avalanches. 
This mutual independence may not be the case for other 
extremal models such as the Sneppen model or Invasion 
Percolation. To understand to what extent the results of 
this work apply to these other models is the direction of 
our current research 

Now we are in a position to write down the exact "mas- 
ter" equation describing how avalanche merging changes 
^(5*, fo) as fo is raised. Let R'^{S, fo) be the average 
number of updated (covered) sites in an /o-avalanche of 
temporal size S. From our simulations of the BS model 
Q we know that for fo close to fc, R'^{S, fo) scales with 
S as S"^^^ , where D is the fractal mass dimension of 
the avalanche. However, for the following arguments any 
form oi R'^{S,fo) will suffice. The "master" equation for 
P{S, fo) can be written as 

(1 - ./o) ^^iy°^ = -p{s, fo)R''{s, fo) 



dfo 



S-l 



P{Si,fo)R''{S,,fo)P{S-Sufo) 



(1) 



Si = l 



Here the first term describes the loss of avalanches of 
size S due to the merging with the subsequent one, while 
the second term describes the gain in P{S, fo) due to 
merging of avalanches of size 5*1 with avalanches of size 
S—Si. It is convenient to change variables from fo to g = 
— ln(l — fo), so that /o = corresponds to g = 0, /o = 1 
corresponds to g = +oo, and dg — This change is 

due to the fact that, although traditionally new random 
numbers are drawn from the flat distribution V{fo) — 
1, the "natural" distribution for the BS model has the 
probability density V{g) = e^^. As usual, the critical 
properties of the model are independent of the particular 
shape of T'. In the rest of the paper we will use the 
"natural" variable g instead of fo- 

To proceed further we make the Laplace transforma- 
tion of P{S,g): p{a,g) = Es=i ^(^'^ 5)e""^ and of 
P{S,g)R''iS,g): r{a,g) = Y.'^^,P{S,g)R%S,g)e~'''' . 
The equation (0) can be conveniently written in terms 
of Laplace transforms p{a,g) = X^sLi ^('5', <?)e^""^ and 
r{a,g) = Es=im5)^'(5,5)e-"^ as dpia,g)/dg = 
-r{a, g) + p{a, g)r{a, g), or simply. 



91n(l -p{a,g)) 



r{a,g) 



(2) 



This exact equation is the central result of this paper. 
It has many interesting physical consequences. When 
g < gc all avalanches are finite (Poo = 0) and nor- 
malization requires p{0,g) = ^(5', .9) ~ 1- From 
the general properties of the Laplace transform one can 



2 



write the Taylor series for p{a,g) and r(a,g) at a = 
as p(a, g) {S)ga + {S^) 30^/2 - {S^)ga^lQ + ... 

and r{a,g) = {R'^)g - {R''S)ga + {R''S^)ga^/2 -+.... 
Substitution of these expressions into Eq. (|^) results in 
^ ln((5),a - (52),aV2 + ...) = {R'')g - {R''S)ga + 
{R'^S^)ga'^/2 + . . .. Since the equation (|) holds for arbi- 
trary a, comparing the coefRcients of different powers of 
a in the above Taylor series results in an infinite series 
of exact equations. Comparison of the coefficients of a" 
gives 



dg 



(3) 



This is exactly the "gamma" -equation derived in 
This equation is valid not only for the BS model but 
for the Sneppen model and Invasion Percolation as well 
since it does not rely on the assumption that the sizes 
of subsequent avalanches are uncorrelated. To show this, 
one has to take a large number N of /o-avalanches and 
write the balance equation of how this number decreases 
as fo is increased |^. Changing the variables from g 
back to fo gives the more familiar form of the "gamma"- 
equation: d\n{S)fJdfo = {R'')fJ{l - fo). 

Higher powers of a in Eq. (||) give new exact equations. 
Here we show just the first two: 

d f{S^) 



dg V {S) 



d {S')g 

dg [ 3{S)g 




2{R''S)g 



{R''S^)g 



(4) 



(5) 



The Taylor expansion changes slightly in the overcrit- 
ical region, where there is a finite probability -Poo (3) 
to start an infinite avalanche. Since the avalanche dis- 
tribution P(5', g) is limited to finite avalanches, it is 
naturally normalized to 1 — Poo{g). Therefore, when 
g > gc the Fourier series for p(a,g) can be written as 
p(a, g) = l-P^ (g) - {S)ga + {S^)ga''/2 + . 
comparison of the coefRcients at in Eq. ( 



dlnPoo(g) 
dg 



.. Now the 
) gives 

(6) 



This new equation is the g > gc analog of the "gamma" - 
equation (||). We will refer to it as "beta" -equation (the 
exponent (5 is traditionally used for the scaling of Poc{g), 
while —7 is used for {S)g.) 

There is a more straightforward way to derive equation 
from the average properties of the merging process. 
The merging of finite and infinite avalanches gives an 
infinite avalanche and, therefore, leads to an increase in 
Pooig). A simple analysis gives an equation governing 
this process as dPoo/dq ~ {R'^)gPoo, which is just Eq. 

®- 

As in the undercritical case, the Taylor expansion of 
Eq.(^ for g > gc determines an infinite series of exact 
equations. Here are the first two: 



dg\PU9) Pcoig)' 



= -{R's), 



= -{R'S^)g 



(7) 
(8) 



As in other creation-annihilation branching processes, 
the avalanche distribution P{S, g) in BS model for g < gc 
is known to have a scaling form 



PiS,g) = S-^fiS'^ig-gc)) 



(9) 



where r and a are some critical exponents and f{x) is 
a scaling function that rapidly decays to zero as x ^ 
—00. From (^) it follows that the average avalanche size 
diverges when g approaches gc from below as {S)g ~ (gc — 
g)~''^ , where 7 = Substitution of this expression into 
the "gamma" -equation O) results in 
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9c - 9 



for g < gc 



(10) 



In the BS model the critical exponent 7 determines the 
amplitude of {R'^)a below gc. The exponent relation de- 
rived from (|lO| ) connects t to D and a: t = 1 + d/D — a 
Q . It is easy to see that Eqs (^jsj) do not yield additional 
exponent relations but further restrict the exact form of 
the avalanche distribution. In fact it can be shown that 
Eq. (H) and the exponents D and t uniquely determine 
the shape of the scaling function f{x) 0. 

The scaling should work in the overcritical regime as 
well. However, unlike in the equilibrium statistical me- 
chanics, the critical exponent a can a priori be different 
above and below the transition. In what follows we show 
that at least for the BS model this is not true. Substitu- 
tion of the scaling form form Poo (5) ^ {9 ^ gc)^ into the 
"beta"- equation (0) results in 



9 -9c 



for g> gc 



(11) 



Again, similar to the "gamma" -equation 



cal exponent f3 gives the amplitude of {R' /g 



the criti- 
above the 

transition. From (11) it follows that the same exponent 
relation t = 1 + d/D — a holds in the overcritical re- 
gion, and, therefore the exponent a is the same above 
and below the transition. The scaling form (^ can now 
be extended to include the overcritical region. For this 
one just lets the argument x of the scaling function f(x) 
be positive. As in various percolation problems j|] the 
scaling form (^) for P{S, g) a,t g > gc results in the ex- 
ponent relation f3 = An interesting consequence of 
exact Eqs. (|l^, [l^) is that the universal amplitude ratio 
for {R'^)g is given by the ratio of two critical exponents 



a-^g 
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1 



2-r 



(12) 
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This unusual relation between the universal amplitude 
ratio and critical exponents is to our knowledge unique 
for the BS model. 

There is a case when the master equation (^) can be 
written in a closed form. This is the extensively stud- 
ied 1^,0 mean field random neighbor version of the BS 
model, where at each time step K — 1 "neighbors" of 
the active site are selected in an annealed random fash- 
ion throughout the whole system. It is easy to see that 
in the thermodynamic limit of this model the number of 
updated sites in the avalanche of temporal duration S 
is given by ricov = [K — 1)5" -I- 1. This is the quantity 
that should be used instead of R'^{S, g) in our equations. 
The missing equation connecting r{a,g) and p{a,g) is 
r{a, g) = —{K — 1) ^^g"'^"* <?) and the final form of 

the Eq. (13) for the mean field BS model is 



d\n{l-p{a,g)) 
dg 



<K-l) 



dp{a,g) 
da 



-p{a,g) (13) 



This equation should be solved with the initial condition 
p(a,0) = e~", since P{S,0) — 5s.i- We checked that 
for if = 2 the generating function Y^^=i P{S^ fo)x 

l-2x/„(l-/„)-[l-4:r/o(l-/o)]^/" 



derived in [TLOl 



ferent methods, after the substitution of x 



using dif- 
e"" and 



e ^ satisfies (113) and has the correct initial 



That confirms the overall consistency of our 



fo = 1 

condition, 
approach. 

It can be shown [Q that the substitution of 
the scaling form (^) into the Eq.(||) defines recur- 
sively all terms in the Taylor series of the scal- 
ing function J[x) at x = 0: /("+i)(0) = 

r(l-T+gni)r(a-+<Tn2) w! f^n^)(r^\f{n■,.)(r\\ 
Z^ni+n2=n r{l-T+cr+an) ni!n2!-' ^ ' 

, where r(a::) is the Euler's gamma function. Actually 
Eq. (H) does more than that: for a given d/D it uniquely 
selects r. We suspect that only for this r the scaling func- 
tion satisfies all usual requirements, such as f{x) —>■ 0, 
when X ±oo, and /q°°x-^(/(0) - fi-x^/''))dx ^ 
(the absence of the infinite avalanche below gc-) Which 
of these constraints actually defines r as a function of 
d/D remains to be determined. The numerical solu- 
tion of the Eq. (0) with R'^iSJo) = AS'^^^ indeed 
seems to give the correct value for r 0|. In Fig. |] we 
present the results of the numerical solution of Eq.(|^) 
with R'^{S,fo) — 5°'^^^, corresponding to the best nu- 
merical estimate of 1 / D in the one-dimensional BS model 
The solution indeed yields r = 1.1 ± 0.1 which is 
consistent with r = 1.07 ± .01 determined by extensive 
Monte-Carlo simulations. 

In 1^ it was shown that all the exponents of a general 
extremal model are determined in terms of just two inde- 
pendent ones, say D and r. As a result of the approach 
described here the set of independent exponents for the 
BS model was narrowed down to just one. It is tempt- 
ing to extend these arguments to other extremal models. 
The weak point for this lies in the assumption that sizes 



of two subsequent avalanches are mutually uncorrelated. 
At present it is unclear how strong is this correlation and 
how it influences the scaling. 
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FIG. 1. The results of the numerical solution of Eq.(l) on 
the interval 1 < 5" < 100 with R'^iS,g) = S°-^^'^. Values of 
g increase from top to bottom. The exponent of power law 
part was measured to be 1.1 ± 0.1. 
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